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DIFFERENCE  METHODS  FOR  STIFF  ORDINARY  DIFFERENTIAL  EQUATIONS 
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Heinz-Otto  Kreiss 


1.  Introduction. 

One  of  the  simplest  examples  of  a stiff  differential  equation  is 
the  scalar  equation 

(i.i)  dy/dt  ^ a, ) y + bedt>  y(°)  = Y0»  t > 0 . 

Here  a^,b,  d are  complex  valued  constants  with  a = Real  a^  ^ < 0 and 
Real  d < 0.  Also  la^  I » 1 but  b/a  and  d are  of  moderate 

4 

size.  A typical  example  is  given  by  a^  ^ = -10  , b/a^  ^ = 1 , d = i. 

The  expression  "of  moderate  size"  is  rather  vague.  It  depends  on  the 
stepsize  k one  wants  to  employ  in  a numerical  calculation.  Often  it 
is  satisfactory  to  say  that  K is  a constant  of  moderate  size  if  Kk  < 0.  l . 
The  solution  of  ( 1 . 1)  is  given  by 

y(t)  = y^t)  + y2(t)  , 

where 


YjU)  = - 


y2(t)  = <y0  + 


Thus  the  solution  consists  of  the  forced  solution  y^t)  and  the  transient 
solution  y2(t).  The  forced  solution  is  smooth  and  varies  slowly.  For 
the  transient  solution  we  have  two  fundamentally  different  situations. 
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1)  a = Reala^  «-l.  In  this  case  y2(t)  decays  rapidly.  Outside 

a boundary  layer  of  order  <&(lal  1 I log  I a 1 I ) we  can  neglect  y^(t). 

1)  Real  a is  of  moderate  size.  Now  y (t)  oscillates  rapidly 
11  ^ 

but  does  not  decay  quickly.  We  shall  not  treat  that  case. 

This  division  of  the  solution  into  a slowly  and  a rapidly  varying 
part  is  typical.  Let  us  consider  a general  system  with  constant  coefficients 
(1.2)  dy/dt  = Ay  + F(t),  y(0)  = yQ,  t > 0, 

where  y = (y(  1 \ . . . , y(n))'  and  F = (F(  1 . . . , F(n))'  1 } are  vector 

functions  with  n components  and  A is  a constant  n x n matrix.  We 
shall  make 

Assumption  1.1.  The  eigenvalues  V of  A can  be  divided  into 
two  sets  M^,  M^. 

1)  X . e M,  if  Real  X « -1  and  | Im  X . I < p I Real  X . I . 

2)  X € M,  if  lx  . I < 7 C . 

j 2 J 2 

Here  p,  C are  constants  of  moderate  size. 

We  need  also  the  following  concepts: 

Definition  1.1.  A matrix  function  A = A(t),  t >0 


A = 


nl 


a 

nn  / 


1 \ 

1 If  y is  a vector  then  y'  denotes  its  transpose  and  y its  adjoint.  The 
vector  norm  is  defined  by  I y | = maxly'1'!-  Similar  notations  hold  for 
matrices,  for  example  I A I = sup  i Ay  (/ 1 y I . 

y 
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is  called  negative  dominant  if  there  is  a constant  p of  moderate  size 
such  that  for  all  t > 0 


1.  3) 


!lm  a | < p I Real  a,. I,  i = 1,2, ...,n  ; 
li  - ii 


and  a constant  6 with  0 < 6 < 1 such  that  for  all  t > 0 

(1.4) 


n 

V 


Rea  a < -1,  a I < - 1 - o)  Real  a , l = 1,  2,  . . . , n . 
n ’ u.  i)  it 

J = 1 

i*J 

It  is  called  essentially  negative  dominant  if  (1.  3),  (1. 4)  are  replaced  by 
(1.  3a) 

(1. 4a) 


Im  a. . < pi  Real  a. . I + c 

ii  - n 


(1.4b) 


Real  a. . < c 
n 

n f-(  1 - M Real  a . + c if  Real  a..  < 0 

\ i i I ii  n 

LJ  13  I -\ 

j = l,  j*i  U t.c  - Real  a.. 


if  Real  a. . > 0 
n 


where  c is  a constant  of  moderate  size. 

It  is  well  known  that  there  is  a transformation  S such  that  the 
matrix  A of  the  system  (1.  2)  can  be  transformed  to 


(1.  5) 


/ A , 0 

A = S'1 AS  - [ 11 

0 A22l 


where  the  eigenvalues  of  A.,  belong  to  M.,  i = 1,2.  A^  is  negative 
dominant  and  |a^^  I < 2c.  Without  restriction  we  can  assume  that  A 
is  already  in  blockdiagonal  form  (1.5),  i.e.  (1.2)  can  be  written  as 

(1.6) 


d_ 

dt 


1 \ 

y 

ii 

\y 

A11  ° 


\ 0 


22 


i \ 

y 

+ 

F 

ii 

Ir-11 

l\y 

1 F 
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is  of  moderate 


Of  course,  we  cannot  guarantee  that  J S I + |S~' 
size.  However,  if  this  expression  is  large  then  the  original  dependent 
variables  y had  not  been  correctly  choosen  and  we  can  consider  S 
as  a scaling  of  y.  In  the  next  section  we  are  going  to  prove 

Theorem  1 . l.  The  solution  of  (1.6)  can  be  written  in  the  form 

y(t)  = ys(t)  + v(t)  • 

Here  yQ(t)  and  its  derivatives  can  be  estimated  by  A~ 1 A . F11 
a 11  22’ 

and  their  derivatives.  v(t)  is  a solution  of  the  homogenous  system  (1.6) 
with  initial  values 

vJ(0)  - y!(0)  - y^o),  vn(0)  = 0 . 

By  assumption  the  eigenvalues  of  A^  have  the  property  that 

Real  « -1.  Therefore  v(t)  decays  rapidly  and,  outside  a boundary 
layer,  the  solution  y of  (1.6)  is  as  smooth  as  yc(t). 

u 

The  last  theorem  is  a special  case  of  the  following  theorem  for  systems 
(1«7)  dy/dt  = A(t)y  + F(t),  t>0, 

with  variable  coefficients.  (See  also  [ 2 ] . ) 

Theorem  1.2.  Consider  the  system  (1.7)  and  assume  that  A is 
essentially  negative  dominant.  Then  the  solutions  of  (1.7)  can  be  written 
in  the  form 

y(t)  = yg(t)  + v(t)  , 

where  yg(t)  and  its  first  p derivatives  can  be  estimated  by  c,  6 , p and 


the  functions 


(i.8)  minfla^l  \ Ijdva  /dtv,  min(|a..|  \ l)d1' F^^/dt1' , v =0,1,2,...,  p. 

v(t)  is  the  solution  of  the  homogenous  equation 
(1-9)  dv/dt  A(t)v,  v(0)  = y(0)  - Yg(0)  , 

which  away  from  a boundary  layer  has  the  same  smoothness  properties 
as  ys(t). 

Assume  that  the  functions  (1.8)  are  of  moderate  size.  Then  the 
boundary  layer  can  always  be  resolved  by  making  a logarithmic  stretching 
of  the  independent  variable  t.  If  we  do  this,  the  solutions  of  (1.7)  will 
be  smooth  everywhere. 

It  should  be  pointed  out,  that  we  do  not  make  any  assumption  that 
the  number  of  "large"  eigenvalues  of  A is  constant.  Thus  we  are  able 
to  treat  tumingpoints . 

One  might  think  that  the  conditions  of  Theorem  1.  2 are  too  restrictive. 
We  shall  give  three  examples  which  show  that  this  is  not  so.  In  all 
these  examples  e > 0 denotes  a small  constant  and  t > -1. 


1)  Consider  the  system 


An  easy  calculation  gives  us 


yU)(t)  = yU)(-l),  t > -1  , 


y(2)(t)  = 


VX)(-1)  + e'(t+i)/e(y(2)(-l)  - y(1)(-D)  for  -1  < t < 0 , 

-L  fe  rl)"rlZdriy(i)(-i)  + e_<j(t’0)y(2)(0)  for  t > 0 , 


where 


ff(t,  r\)  - e 


J e“^'/ede,  i.e.  0 < <r(t,  r,)  < 1 . 


Thus  y(2\t)  changes  rapidly  in  a neighbourhood  of  t-0  and  becomes 
of  order  6(e  'i/  I y^  ^(-1)  I)  f°r  t > 0. 


2)  Consider  the  system 


0 e 


a^lt)  e 


. y,  y(~f)  = * 

A r\ 


where 


0 for  t < 0 , 


-t  for  t > 0 . 
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Now  y(t)  = (1,0)'  for  t < 0 and  converges  rapidly  to  zero  for  t >&(\e). 

Observe  that  the  coefficients  a or  a *da . ,/dt  are  bounded. 

ij  U iJ 

Thus  the  diagonal  dominance  of  A is  important. 

3)  One  might  think  that  one  has  only  to  control  the  eigenvalues  of 


A and  their  variation  to  be  sure  that  the  solution  of  the  differential 
equations  is  of  moderate  size  and  smooth.  This  is  not  so.  Consider 


the  system 


dy  * * n \ 

e T~  - U (t)  U(t)y  = A(t)y,  y(0)  = y , 

dt  0-1  ° 


where  n > 0 is  a constant  and 


cos  at  sin  at \ 

U(t)  - , a = const. , 

-sin  at  cos  at 


is  a unitary  transformation.  Let  v = Uy  denote  new  dependent  variables. 
Then  v is  the  solution  of 


-1/e  n/e  - a^ 


v = Bv  . 


This  is  a system  with  constant  coefficients  and  its  general  solution 
has  the  form 

KjX  K^X 

v(x)  = e / + e k . — — e ± Na(n/e  - a)  . 

l z ) 

Here  are  the  eigenvalues  and  t.  are  the  corresponding  eigenvectors 

-»  -2 

of  B.  As  long  as  a(g/e  - a)  < e the  solutions  of  the  system  decay 


exponentially.  However,  if  a (n/e  - a)  > e then  there  is  one 


-7- 


. 


exponentially  increasing  and  one  exponentially  decreasing  solution. 

-[ 

This  happens,  for  example,  if  a - 2,  n e >2.  Observe,  that 
the  eigenvalues  of  A(t)  are  constant  and  that  U(t)  is  a slow  rotation. 

The  above  examples  show  that  one  needs  to  be  cautious  when 
solving  stiff  systems.  It  is,  of  course,  not  necessary  that  the  system 
(1.7)  be  essentially  diagonally  dominant.  It  is  sufficient  that  one  can 
transform  it  to  this  normal  form.  In  section  three  we  shall  show  that 
this  is  always  possible  using  an  adequate  stretching  of  the  independent 
variables . 

We  shall  now  discuss  difference  approximations.  We  divide  the 
t-axis  into  subintervals  of  length  k,  define  gridpoints  t^  = vk, 
v - 0.  1.  2,  . . . and  denote  by  u = u(t  ) vector  functions  defined  on 

777  V V 

the  grid.  We  approximate  (1.7)  by  a muitistep-method  which  is  generated 
from  an  identity 

r 

(1.10)  £ aJy(tv.j)+  kpjdy(ti)_.)/dt  - kv|»(tv)  , 

where 

(1.11)  440  = <$(kr  1dry/dtr) 

denotes  the  truncation  error.  Observe  that  (1.10)  has  nothing  to  do  with 
the  differential  equation.  It  is  valid  for  all  sufficiently  smooth  functions. 

Only  when  we  replace  dy/dt  by  Ay  + F the  differential  equation  enters 
and  we  obtain  the  corresponding  multistep-method, 

i 
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exists  if  the  eigenvalues  K of  A satisfy  the  inequality 

kp  ^ Real  X > -1  • 

Now  the  solution  of  (1.12)  is  uniquely  determined  if  we  specify  initial  values 

u - y p=0,l,2,...,r. 

p p M-  P 

The  approximation  is  only  useful  if  it  is  stable.  However,  here  we 
cannot  define  stability  by  describing  the  behavior  of  a method  as  k - 0, 
because  we  want  to  use  it  in  the  case  that  k j A j » 1.  Instead  we 
consider  classes  of  problems  (1.7),  and  define  uniform  stability. 

Let  £2  be  a domain  in  the  complex  plane  with  the  following  properties 

1)  There  is  a constant  C (of  moderate  size)  such  that  z « £2 
and  Real  z >0  implies  I z | < C . 

2)  If  z ( (2  then  also  jz  e £2  for  all  real  a with  0 < a <1. 

We  shall  consider  the  class  ft  (£2)  defined  by 

Definition  1.  2.  7r(£2)  = ^(£2,  c,  p,  6,  K)  denotes  the  class  of  problems 
(1.7)  where  A is  essentially  negative  dominant,  the  functions  (1.8)  are 
bounded  by  a constant  K (of  moderate  size)  for  p = 1 and  the  eigenvalues 
X of  A belong  to  £2. 


-9- 


Stability  is  defined  almost  as  usual.  We  assume  that  one  can  solve 


(1-12)  for  u boundedly  and  that  the  solutions  of  (1.12)  increase  at 
most  slowly.  The  only  difference  is  that  all  the  constants  involved  are 
independent  of  the  particular  problem  but  hold  for  the  whole  class. 

Definition  1.  3.  Consider  a class  v of  problems  (1.7).  A multi- 
step-method (1.12)  is  stable  for  the  class  if  there  are  constants  K , K , 

J b 

, (0) 

and  k >0  such  that  for  all  problems  in  y , all  t,  and  all  k 
writh  0 < k < k^ 

(1-13)  |(a_1I  + kp_1A(t))"1|  < , 

and  the  solution  of  the  homogenous  equation 

r 

(1.14)  L[vJ  = V (a. I +k(3.A(t  ))v  . = 0 

j__1  1 J v-}  v-) 

satisfy  the  estimate 


(1.15) 


\c(t  -t  ) r 
Iv(t  )l  <K„e  S * v X 

11  s j = 0 


t > t . 

M-  V 


Already  G.  Dahlquist  [1]  has  defined  the  stability  for  classes  of 
problems.  He  considered  the  class  7 ? of  all  scalar  equations 

(1.16)  dy/dt  = ky,  \ = const.  Real  k <0  , 

and  proved  that  the  trapozoidal  rule  is  the  "best  method"  (method  of  the  highest 
order  with  the  smallest  coefficient  of  the  truncation  error)  which 
is  stable  for  the  class  7?  . Later  O.  Widlund  [8]  introduced  the  class 
7?  = ?/j(p)  of  all  scalar  equations  (1.16)  where  k satisfies  the  condition 


-10- 


tl-17) 


Real  X < 0,  |lm  X i < p I Real  X- 1 , 


and  showed  that  there  are  methods  of  hiyher  order  than  two  which  are 
stable  for  such  a class.  There  are  now  a large  number  of  stable  methods 
available.  See  for  example  l 5 ] , [9].  More  generally,  one  can 
consider  classes  r(Q)  of  scalar  equations  (1. 16)  where  X belongs 
to  a domain  Id  of  the  above  type. 

We  consider  now  the  approximation  (1.12)  for  the  class  . In 

this  case  the  solutions  can  be  given  explicitly.  They  are  of  the  form 

(1.18)  uv  = ^ Pn(v)V 

where  the  k are  the  solutions  of  the  characteristic  equation 
P 


(1.19) 


r 

^ (a . + Xkfl  ,)k  1 = 0 

. LJ,  ) ) 

; = -l 


and  the  P are  polynomials  in  v of  order  one  less  than  the  multiplicity 
P 

of  k . We  make  the  following 
P 

Assumption  1-2.  1)  The  approximation  (1.12)  is  stable  in  a neighbour- 

hood of  Xk  = 0 which  means: 

a)  For  Xk  = 0 the  solutions  of  (1.19)  satisfy  the  condition 

(1.20)  I k (0)  I £ 1.  If  U (0)|  = 1 then  its  multiplicity  is  one  . 

P P 

b)  No  weak  instability  occurs,  i.e.  in  a neighbourhood  of  Xk  = 0 

the  solutions  k (Xk)  with  Ik  (0)  | = 1 are  of  the  form 
P P 


(1.  21) 


k (Xk)  = k (0)(1  + g Xk)  + 0(Xk)  with  g > 0 . 
u p P R 


-11- 


A 


1)  The  approximation  is  stable  in  a neighbourhood  of  \.k  - 00  . 

(Dividing  ^1.19)  by  kk  we  can  consider  its  solutions  as  functions  of 

(\k)  * and  the  stability  is  defined  in  the  same  way  as  above.) 

3)  Let  d ,d.  with  0 < d < d < °°  be  constants  and  denote  by 

Q , , the  domain 

dl’d2 

z e U , , if  Real  z < 0 and  d < |z  I £ d . 

drd2  ~ 1 d 

If  \k  « n fl  Q , , then  the  solutions  k of  (1.19)  satisfy 

dl’d2 


(1-  22) 


U | < 1 - t,  t > 0 a constant  depending  on  d , d . 
P 1 


Remark.  The  last  condition  is  often  strengthened  to: 

For  every  fixed  d^  > 0 there  is  a constant  t such  that  (1.  22) 
holds  for  all  In  this  case  we  call  the  method  strongly  stable. 

As  an  example  we  consider  methods  which  are  generated  from  the 
identity 

(1.  23)  y(t  ^ - akdy(t  ,.)/dt  = y(t  ) + (1  - a)dy(t  )/dt  + M>(t  ) , 

v tI  w+I  v v v 

where  a is  a constant.  (1.2  3)  leads  to  the  approximations 

1 + k(l  - o~ ) k 

Uv+1  = 1 - ka\  Uv  ’ 

which  are  unstable  for  any  class  7/j(p)  if  <r  stable,  but  not 

strongly  stable,  if  <x  = ~ and  strongly  stable  if  <r  > j . 

Now  consider  a class  Y of  problems 


dy/dt  = A(t)y  + F 


-12- 


and  assume  that  the  class  is  such  that  the  eigenvalues  \(t)  of  every  A 
belong  to  a set  U.  Approximate  these  problems  by  a multistep-method 
which  is  stable  for  the  class  7?(Q).  The  main  q"estion  which  we  want 
to  discuss  in  Section  4 is  whether  the  multistep-method  is  also  stable 
for  the  class  v.  In  this  generality  the  answer  is  no,  as  we  shall 
demonstrate  now. 

Consider  our  third  example,  i.e. 

* /_1  n\ 


(1.  24) 


• £ ■ 


U(t)y  - A(t)y  , 


\0  -1/ 

and  approximate  it  by  the  multistep-method  generated  by  (1.  23)  with 
1 

o-  > g . 1 • e • 


(1.  25)  (I  - kcrA(tw+1»Uv+1  = (I  + k(l  - <r)A(t^))uv  . 

Introduce  into  (1.25)  new  variables  w = Uu.  Then  we  obtain  a system 


with  constant  coefficients 
(1.  26) 

where 


w , = Bw 

v+1  V 


B = I-2^ 

£ 


f 1 r|  cr 

10  1 


-1  r) \\  / cos  ak^  sin  ak^ 


0 -1 1/  \-sin  ttk^  cos  ak^ 


I + (1  - cr)k 


-1  n 

0 -1 


o 

l\ 

\l 

I + ak 

1 

1-1 

o, 

+ G(ak)  1 

p n(l  - o-) 
o P 


L 


(1.  26)  is  an  approximation  of  the  homogenous  system 


(1.  11) 


d v 
dt 


v = Bv 


-1  r)  - ae^ 
at  -1 

and  we  know  that  the  solutions  of  (1.  27)  decay  exponentially  as  long  as 
ear)  = t < 1.  We  shall  now  consider  three  cases. 

1)  cr  = -j,  e « k.  A somewhat  tedious  calculation  shows  that 
eigenvalues  k of  B are  approximatively  given  by 


r 2~z 

V 8 a q£  -a  k 


If  a > 0 and  e = ~ ak^  then  we  get  for  the  negative  root 
4 


k » —1  + ak(l  - n/2 r|  - 1)  , 

and  |k  I > 1 for  r|  > 1.  Thus  the  approximation  has  an  exponential 

increasing  solution,  though  the  solution  of  the  differential  equation  decays 

2 2 

as  long  as  q < 1/ea  = 4/a  k . Furthermore  the  increasing  exponential 
solution  can  grow  arbitrary  fast  and  therefore  the  method  is  not  stable 
for  the  class  of  problems  of  type  (1.24)-. 

2)  <r  = 1,  e « k.  Now  we  get 

t f[  2 

k « - “ ± S~  t t = ear|  . 

2 4 

If  |t  | < 1 the  method  is  stable.  However,  if  t < -1,  then  the  method 
is  not  stable  though  the  solutions  of  the  differential  equations  decay 
exponentially  as  t — ».  (q  = ©(e  *)  and  the  matrix  B is  far  from 
being  negative  dominant. ) 
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3)  ■“  < o < 1,  e « k.  In  this  case  the  approximation  has  an 
exponentially  increasing  solution  lor  lean  I > tq  where  is  some 

constant  with  0 < t < 1.  Again  the  matrix  B is  far  from  being 
negative  dominant. 

In  some  sense  the  most  disturbing  result  is  the  effect  on  the  trapezoidal 
rule  (cr  - ) because  the  difference  approximation  has  exponential 

increasing  solutions  for  examples  which  cannot  be  considered  pathological 
at  all.  This  kind  of  behavior  is  typical  for  methods  which  are  only  stable 
and  not  strongly  stable  for  a class  7?(Q). 

The  last  example  shows  that  one  cannot  decide  the  stability  of  a 
method  for  a given  class  of  problems  by  just  looking  at  scalar  equations. 

One  has  to  impose  other  conditions.  The  main  result  of  Section  4 is 
Theorem  1.3.  Consider  a class  fr(n)  and  assume  that  the 
approximation  (1.12)  satisfies  assumption  1.2  for  the  corresponding  class 
??(q).  Then  it  is  stable  for  the  class 

We  shall  now  derive  error  estimates.  Consider  the  problem  (1.7) 
and  assume  that  it  has  a smooth  solution  y(t),  i.e.  a number  of  its 
derivatives  are  of  moderate  size.  If  the  conditions  of  theorem  1.  2 are 
fulfilled  and  the  functions  (1.8)  are  of  moderate  size,  then  this  is  no 
real  restriction.  Either  we  change  the  initial  conditions  to  yg(0)  or 
we  perform  a logarithmic  scaling  of  t which  resolves  the  boundary  layer 
so  that  also  the  solution  of  (1.9)  is  smooth.  We  approximate  the  problem 
(1.7)  by  the  multistep-method  (1.12).  Introducing  y(t)  into  (1.12)  gives 


us,  using  (1. 10) 


L[  y 1 = kG  + k4<(t  ) 

V V 


Therefore  we  get  for  the  error  w = y - u 


(1-^6) 


L[  w ] - k_4>( t ) w = <fi  = 0(k  ),  p.  = 0, 1,  . . . r 

v M-  M- 


If  the  approximation  is  stable,  then  (1.26)  gives  us  the  usual  error 
estimate.  Observe  that  ij>  depends  only  on  the  derivatives  of  y(t) 
and  not  on  the  coefficients  of  the  differential  equation.  Thus  the  stiffness 
does  not  enter.  Thus  we  need  only  to  investigate  the  smoothness  of 
the  solution  of  the  differential  equations  and  to  derive  stability  conditions. 

Instead  of  starting  from  the  relation  (1.10)  we  can  start  from 
identities  which  contain  higher  derivatives,  for  example 


(1.27) 


y(tv+i)  - ikdy(tv+i)/dt  + II  ^d2y(tv+1)/dt2  = 


y(t  ) + 7 kdy(t  )/ dt  + ~ k2d2y(t  )/dt2  + k^Jj  (t  ) , 
v 2 v 12  v 5 v 


(1.28) 


- f kdy(t„+i,/dt  + i k'd2y<tl,+1)/dt2  = 

y(t„)  + 3 kdylO/dt  + k^ltj  , 


where 


U5(tv/)  I < ~ max  |d5y/dt5|,  U4(tv)l<^  max  |d4y/dt 


1 <t<t 
v — v -fl 


V v+1 
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Observe  the  extremely  small  constant  in  front  of  d y/dt  . These  two 
examples  are  special  cases  of  the  Fade' approximations  (p  = q = 2,  and 
p 2,  q ; 1 respectively) 

\ (-l)Jr!(r  t g - J ! . J w,J  _ 


(1.  29) 


\ rMr  t 4 j j , )/d  > 

J Mr  - j) ! (r  + q)!  U ^ 

1 = 0 


■?  qnr  * g - 1)! k)d)v,.  ,/dt)  + kr+q+14;  (t  ) 

A ,!(r+q)!(q-j)!  kdy(‘w)/dt  Vqtl'V 


where 


K+q+l(t^  * - (r  + q)!(r  + q + 1)! 


I r+q+1  r+q+1 

max  |d  y/dt 

t <t<t 

V V +1 


If  y(t)  is  the  solution  of  the  differential  equation  (1.7)  then  we 
can  rewrite  the  above  relations  as  relations  between  y(t^+j)  and  y(t^). 
The  coefficients  then  depend  on  A,  F and  their  derivatives.  We  obtain 
the  difference  approximations  by  replacing  y(t^+1)  and  y(t^)  in  these 
last  relations  by  u , and  u , respectively,  and  neglecting  4*. 

For  example,  if  y(t)  is  the  solution  of  the  scalar  equation  (1.21)  then 

dJy/dtJ  = ^y 

and  the  difference  approximations  corresponding  to  (1.27)  and  (1.28)  are 

U - 2 kK  + II  ^K+l  = (1  + 2 kK  + 12  kV>Uv 


(1  + jk^uv+i  = ^ " 3 kk  + 6 kV)Uv 


respectively. 
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Lhle  [ i|  has  proved  that  the  methods  based  on  (1.  29)  are  stable 


for  the  class  if  p = q and  strongly  stable  if  p = q + 1 or 

p - q +2.  There  are  no  new  difficulties  in  proving  theorem  1.  3 also 
for  this  type  of  approximations. 

There  is  a large  literature  on  methods  of  this  type,  for  example 
l ^ ] , [ 7 ] , [ 9 ] . Sometimes  they  are  not  derived  from  identities  between 
smooth  functions  and  their  derivatives.  This  can  be  dangerous. 
Approximate,  for  example,  the  differential  equation  (1.1)  by 

U-  30)  (1  - ka  )4u  = (1  - 3ka  )u  + k(l  - ka  )4bedt  , 

11  y+1  11  v 11  ’ 

which  is  strongly  stable  for  the  class  7?  . It  is  consistent  in  the  usual 

sense,  i.e.  the  solutions  of  (1.  30)  converge  to  the  solution  of  the 

differential  equation  (1.1)  as  k -►  0.  However,  the  convergence  to  the 

smooth  part,  y (t)  is  not  uniform  for  the  whole  class.  Let  b = a,, 

11 

and  yQ  + b/(a^  - d)  = 0.  Then 

d t 

lim  y = -e  but  lim  |u  | = oo  . 

aii^-°°  air-°° 

There  are  no  difficulties  to  derive  algebraic  conditions  for  uniform 
convergence.  One  can  develop  the  solutions  of  general  systems  (1.7) 
into  asymptotic  series.  The  same  is  true  for  the  solutions  of  the  difference 
approximations . Comparing  these  series  gives  algebraic  conditions  for 
uniform  convergence. 
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2.  Thu  analytic  problem. 

In  this  section  we  consider  the  initial  value  problem  (1.  7) 

(2. 1)  dy/dt  - A(t)y  + F(t),  y(0)  = yQ,  t > 0 , 

where  Ait)  is  negative  dominant.  This  is  no  restriction,  because  if 

A(t)  is  only  essentially  negative  dominant  then  we  introduce  a new  variable 


~ -at 
y - e y 


and  obtain  the  system 


~ ~ -at 

dy/dt  = (A(t)  - al)y  + e F , 


which  is  negative  dominant,  provided  a is  sufficiently  large.  We  want 
to  estimate  the  solutions  of  (2.1)  and  start  with 

Lemma  2.1.  Consider  the  differential  equation 
(2.2)  dy/dt  = an(t)y  + F(t),  a(t)  = Real  a^  < 0 . 


Then  the  following  estimates  hold 


(2.3) 


(2.4) 


| y(t)  I < t max  I F(  rj)  I + s(t)  I y(0)  I , 

0 <n  <t 

I y(t)  I < max  )F(  rT)/a(n>  j + S(t)  I y( 0)1  if  a < 0 , 

o<n<t 


where 


s(t)  = e 


/ a(£)d£ 
0 


Proof.  The  solution  of  (2.  2)  can  be  written  down  explicitly 


y(t)  = / 

0 


t / au(g)d| 


f a^£)d£ 

F(rj)dn  + e°  y(o)  . 
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The  first  inequality  follows  from 


t t 

/ au($)d!  / a(£)d£ 


< e 


< 1 . 


Furthermore  if  a < 0 then 

. /v6,ds 

i y(t)  I - s(t)  I y(0)  | <1/  e n F( rj)d -q  I £ 

0 

t 

t / a(£)d£ 

I I a(n)e  n ^(inJ/aln)  i d Tq  | < 

0 

t 

t / a(t)d£ 

I f ~ e n d n I ' max  i F( n)/a(  n)  I , 

o dn  o < n <t 

which  gives  us  (2.4). 

We  consider  now  the  homogenous  system 
(2.5)  dv/dt  = A(t)v  . 

For  its  solutions  we  have 

Lemma  2.2.  Let  tQ,  t be  real  numbers  with  0 < tQ  < t.  Assume 
that  A is  negative  dominant  then 

t 

6 / a(|)d£ 

|v(t)|<e  ° lv(t  ) I,  a = max  Real  a..  , 

i 

i.e.,  the  solution  operator  S(t,tg)  of  (2.5)  satisfies  the  estimate 
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(2.6) 


6 / a(t)dt 


I S(t,  tQ)  1 < e 

Proof.  Let  k = (t  - t )/N,  N natural  number,  and  denote  by 
t + vk  gndpoints  and  by  = w(t^)  functions  defined  on  the  grid 


l0  h 


tN-l  tN 


We  approximate  (2.5)  by 


(I  - kAltv+l))Ww+l  = w,'  w0  = v,t0)  ■ 

By  assumption  A is  negative  dominant.  Therefore  for  sufficiently  small  k 

(1-  6ka(tv+1))|wv+1l  < IwJ 


IN-1  i 

< TT  (1  - 6ka(t^+1))  1 WQ 


We  know  that,  as  k -*  0, 


6 / a(|)d£ 


N-l  _L  t0 

w^  - v(t),  TT  (1  - 6ka<tl/+1))  e 

v - 0 


Therefore  the  estimate  (2.6)  follows. 
For  the  system  (2.1)  we  have 
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4 


Lemma  2.3.  Assume  that  A is  negative  dominant.  Then  the 
solution  of  (2.1)  satisfies  the  estimates 
(2-7)  ly(t)  I < t max  | F( g)  | + s( t)  | y( 0)  | , 

0 < T)  <t 

<^.a)  I y(  t)  1 < 26  1 max  !( A(  n)  + A*  r1))"1F(  n)  J + s(t)  | y(0)  | , 

°<n<t 

where 


and 


t 

6 f a(£)d£ 
s(t)  = e ° 


max  Real 
i 


a. . 

it 


9 


Proof. 


hi 

0 . 

. . 

■ • 0 \ 

A = 

0 

a22 

0 . . 

. . 0 

U . 

• 

. 

nn 

The  first  estimate  follows  from  lemma  2.  2 and  the  representa- 


tion 


t 

y(t)  = f su.nJF^Jdn  + S(t,0)y(0)  . 

0 

We  are  now  going  to  prove  the  second  estimate.  By  lemma  2.  2 we  can 
assume  that  y(0)  = 0.  Let  t be  fixed  and  denote  by  M the  space 
of  all  continuous  functions  g(t)  with  g(0)  = 0.  The  system 

S-qY  = dy/dx  - Ay  = F 
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has  a unique  solution  in  M and,  by  lemma  2.1 


12.9) 


1 r01  F il  < 2 lit  A + A*)  1 F || , ilg 


max  ig(£) 

o < i <t 


Let  S*^  denote  the  operator 


S y - (A  - A)y,  y f M . 

Then  we  can  write  the  differential  equation  (2.1)  in  the  form 

(I  - f'^y  = S.QlF  . 

Definition  1.1  and  (2.9)  imply 

Hx'^y  ||  < 2 j|  ( A + A ")  ^y  II  £ (1  - &)  II  y II  , 

and  the  estimate  (2.8)  follows  from  (2.9). 

The  assumption  that  A is  negative  dominant  implies  that 

•1.  : . - 'A-1*  ! ■ I'  * A *'_1(A  - A)  I < ^ (p  + 1 - 6) 


2.10)  |(A  + A ) A I < ! ( A • A ) A I + i(A  + A ) 


Thus  (2.8)  can  also  be  expressed  as 


(2.11) 


ly(t)  I < 


+ 1-6 


|a  \n)F(r|)  I + s(t)  ly(o)  I . 


max 


,[ 


o < n <t 

We  shall  now  prove  theorem  1.2.  We  use  the  notation  F 1 = d F/dt 
and  start  with 

-1  A[  v J 

Lemma  2.  4.  Assume  that  A is  negative  dominant.  A A , 

^-lpt  v 1 , i/  = o,  1,  2,  . . . , p;  can  be  estimated  by  the  functions  (1.8). 

Proof.  Denote  by  A the  diagonal  of  A.  Then 

|A_1A[  | < 1(1  + A_1(A  - A))-1 1 |A_1A[  ^ I • 
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By  assumption  la.  J >1  and  1(1  + A ‘ A - A))  i < 6 Thus  A A 

can  be  estimated  by  A *A^  ‘ which  is  composed  of  terms  a ^ . 

it  ij 


The  same  procedure  works  for  A 


-1 J v 


and  the  lemma  is  proved. 


Lemma  2.  5.  Assume  that  A is  negative  dominant  and  define  C(t) 


by  Aft)  A(0)C(t).  Then  Cl  (C  )' 


0, J,  2,  . . . , p;  can  be 


estimated  by  A A1  , u-°,f>2,...,p. 

— 

Proof.  Let  B A dA/dt.  Then  A is  the  solution  of 

dA  '/dt  = bV  . 

C (t)  is  the  solution  operator  of  this  differential  equation  and  therefore 
the  lemma  follows  from  well  known  estimates. 

Lemma  2.  6.  Assume  that  A is  negative  dominant.  The  solution 
y(t)  of  (2.1)  and  its  first  p derivatives  can  be  estimated  by 


(2.12) 


y[v](0),  A~lA[v],  A_1F[w],  v = 0,  1,  2,  . . . , p 


Proof.  Differentiating  (2.1)  gives  us 


(2.13) 


dyL  '/dt  = Ay1  1 + 


M . v x 


y + FL 


Therefore  the  lemma  follows  by  induction  from  (2.11). 

A simple  consequence  of  lemma  2.  6 is 

Lemma  2,7.  Assume  that  A is  negative  dominant  and  that 


2.14) 


y(0)  = 0,  F (0)  = 0,  v = 0,1,  2,.  . .,p  - 1 . 


Then  the  solution  y(t)  of  (2.1)  and  its  first  p derivatives  can  be 
estimated  by 


(2.  15) 


v 0,  1,  2,  . . . , p . 


A-iAM;  a-1f[v]> 

Proof.  (2. 13)  implies  y^(0)  = 0 for  i/  = 0, 1,  2,  . . . , p.  Therefore 
the  lemma  follows  from  lemma  2.6. 

We  can  now  prove  the  first  part  of  theorem  1.  2 by  reducing  the 
general  case  to  the  special  case  of  lemma  2.  7.  Consider  the  system 
i2.16)  dw/dt  = Aw  + F , 

where 

~ P~1  ^ i i „ P-1  v , , 

A = >,  ~ A V (0),  F = X — FM(0)  , 

v'  = 0 v - 0 

in  a neighbourhood  of  t - 0.  We  seek  a solution  of  the  form 

p-1  V 

(2.17)  w(t)  = ^ ^77  Ww  + <p(0)  = 0 . 

y - 0 


Introducing  (2.17)  into  (2.16)  gives  us 

(2.18)  d<p/dt  = A<p  + f(t),  <p( 0)  = 0 , 

with 


f(t)  - 


p-i  ^ 

L-J  v I 
12  = 0 • 


i Mo> 

■ u~  o 


v ! w 


1 *i_"  - w + F ^(0) 


p!(  v - p) ! v+1 


+ , 


where  w^  = 0 and  f^(t)  is  a polynomial  in  t with  coefficients  depend- 
ing linearly  on  A^1  ^(0)w  . Choose  now  the  w such  that 

P V 


v , . v fw 

(2.19)  y.  A ^ (0) ^ 


! ( v - p) ! v+1 


- w . +F[V)(0)  = 0,  v = 0,1,2,  ...,p  - 1 . 


We  shall  show  that  this  is  always  possible.  Let 


V . r i v ! w 

V A \o)A^'  


i = l 


.\(v  - p)f 
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Then  (2.19)  can  be  written  in  the  form 


U-20)  w„  - A'^OJw^j  = L^  - A_1(0)rlv,(0),  v = 0,l,2 1;  LQ  = 0. 

Observe  that  L depends  only  on  w with  u < v and  that  w =0. 

l'  P P 

Therefore,  if  we  neglect  the  term  A \o)w  , then  the  w are  uniquely 

v +1  V 

ietermined  by  (2. 20)  and  can  be  estimated  by  A *(0)A^  *^(0),  A \o)F^  ^(O). 
The  same  is  true  if  |A  (0)1  is  sufficiently  small.  If  A *(0)  is  not 
small,  then  the  determinant  D of  the  linear  system  (2.20)  can  vanish. 
However,  by  introducing  a new  variable  y = exp(at)y  into  (2.1),  we 
can  always  choose  a such  that  D#0.  Observe  that  D does  not 
vanish  identically. 

In  a neighbourhood  of  t = 0 the  system  (2.18)  satisfies  the 
conditions  of  lemma  2.7.  Furthermore  the  derivatives  A-iA^,  A~M  ^ 
can  be  estimated  by  A ^ V ^ , A M v ^ . Therefore,  also  <p( t)  can  be 
estimated  by  these  functions.  This  is  also  true  for  the  solution  of  (2.16) 

if  we  choose 


(2.21) 


w(0)  = w , i-  e. , <p( 0)  = 0 


as  the  initial  value. 


Let  g(t)  « C be  a "cut-off"  function,  i.e.  there  are  constants 
£*r  ° £ with  0 < « < « < - such  that 

S’ 

1 for  0 < t < , 


g(t)  = 

Let  w be  the  solution  of  (2.16).  Then  w = gw  is  the  solution  of 


0 for  t > a 2 • 
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U-  22) 


dw/dt  = Aw  + F,  F (dg/dt)w  + gF,  w(0)  = 


If  we  choose  the  a j - 1,  2;  sufficiently  small  then  w has  the  same 
properties  as  y in  lemma  2.7.  Subtract  (2.  22)  from  (2.1).  Then 
u - y - w is  the  solution  of 


;-2- 2 3) 


du/dt  - Au  + (A  - A)w  + F - F,  u(0)  = y(0)  - w^  . 


u can  be  written  as  u = u^  + v where 

dv/dt  = Av,  v(0)  = y(0)  - w , 

and  Uj  is  the  solution  of  the  differential  equation  (2.23)  with  homogenous 
initial  value.  For  the  conditions  of  lemma  2.7  are  fulfilled  and 
the  first  part  of  the  lemma  is  proved. 

To  prove  the  second  part  of  theorem  1.2  we  need  another  version 
of  lemma  2.  7. 

Lemma  2.7a.  Assume  that  A is  negative  dominant  and  that 

y(°)  = o,  F^(0)  = o,  V = 0,  1,2,  . . . , p - 2 . 

Then  the  solution  of  (2.1)  and  its  first  p derivatives  can  be  estimated  by 

aV"*,  a'M1'*,  i/  = 0,1,  2,  . . . ,p;  and  F^p-1^(0). 

Proof..  (2. 13)  implies  y^  ^(0)  = 0,  v = 0, 1,  2,  . . . , p;  and 
y^(0)  = F^  P 1(0).  Therefore  the  lemma  follows  from  lemma  2.6. 

We  need  also 


Lemma  2.8.  Assume  that  A(t)  is  negative  dominant.  Let  v(t)  be  a 
solution  of  (2.  5).  Then 
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!tlv)l(i  , t"v  v >p,  = 0,1, 2, ...,p  ; 


can  be  estimated  by 


Let 


A t _ 0, 1,  2,  . . . ,p;  and  v(0)  . 

Proof.  It  is  clear,  that  we  need  to  prove  the  theorem  only  for  ^ 

V 

w - t v.  Then  w is  the  solution  of 


dw/dt  = Aw  + wtl  *v,  w(0)  = 0 . 

We  have  to  prove  that  w^  ^ ^ , p = 0, 1,  2, . . . , v;  can  be  estimated  by 
-1  [ T ) 

A A , t _ o,  1,  2,  . . . , p.  For  p = 1 this  follows  directly  from  lemma  2. 
because  A dv/dt  - v.  The  general  case  follows  by  induction  using 
lemma  2.  6,  and  the  observation  that  v^  ^ can  be  expressed  by  deriva- 
tives of  lower  order  using  the  differential  equation  (2.  5). 

The  last  lemma  shows,  that  outside  a "small"  neighbourhood  of 

t = 0,  the  solution  v(t)  of  (2.  5)  and  its  derivatives  can  be  estimated 

-1  [ v] 

by  A A . We  shall  now  improve  these  estimates.  We  shall  show, 
that  if  A dA/dt  is  of  moderate  size,  then  the  rapidly  changing  part  of 
v(t)  is  to  a first  approximation  the  solution  of  the  system 
(2.24)  dv^/dt  = A(0)Vj,  v^(0)  = v(0)  , 

with  constant  coefficients. 

Theorem  2.2.  Assume  that  A is  negative  dominant.  Then  the 
solution  of  (2.5)  can  be  written  as 
(2.2  5)  v = vi  + ei 
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where  v^  is  the  solution  of  (2.24)  and  e^  is  the  solution  of 
(2.26)  deA/dt  - A(t)e1  + F(t),  F(t)  = (A(t)  - A(0))  v^t),  e^O)  = 0 . 

furthermore  e^  and  de^/dt  can  be  estimated  by  A ^dA/dt. 

Proof.  Introducing  (2. 25)  into  (2.  5)  gives  us  (2.  26).  By  lemma  2.  5 

A~\t)F(t)  = 1 ~ ^ A-1(0)tv  ■ 

Lemma  2.8,  applied  to  the  equation  (2.  24),  and  lemma  2.5  show,  that 
the  theorem  follows  from  lemma  2.  6 with  p = 1. 


In  the  same  way  one  can  prove. 

Theorem  2.3.  Assume  that  A is  negative  dominant.  Then  the 
solution  of  (2.5)  can  be  written  as 

XP 

(2.  27)  v(t)  = 2,  v (0  + e (t) 

P P 

where  v (t)  is  the  solution  of  (2.  24),  the  v (t),  pi  > 1 are  the  solutions  of 
1 P 

dv  Vdt  = A(0)v  ,.(t)  + (Aft)  - A(0))v  , v (0)  = 0 , 

p.+l  ^i+l  ^ 

(2.28) 

p = 1,  2, . . . , p - 1; 

and  e satisfies  the  equation 
P 

(2.29)  dep/dt  = A(t)e  + (A(t)  - A(0))vp(t),  ep(0)  = 0 . 


Furthermore 


0, 1,  2,  . . . , p, 


0,1,2,  . . . , p. ; p.  - 1,2,  . . . , p - 1 , 


-29- 


d 


can  be  estimated  by 

v(0)  and  A *A^  ^ , jx  = 1,  2,  . . . , p . 

Finally,  in  the  same  way  as  before,  we  can  replace  A(t)  by 
P-1  tv  r | 

^ ~ A 1 (0)  without  destroying  the  desired  estimates.  The  resulting 

v - 0 v ' 

equations  can  be  solved  explicitly  and  consist  of  terms  of  the  form  ^e^a 
where  denote  the  eigenvalues  of  A(0).  These  terms  generate  at 
most  boundary  layers.  We  have  thus  proved  theorem  1.2. 


r 

i 


3.  A normal  form  for  the  differential  equations. 

In  the  first  section  we  discussed  examples  which  show  how  important 
it  is,  that  the  system  (1.7)  is  essentially  negative  dominant  and  that 
the  functions  (1.8)  are  of  moderate  size.  One  can  also  give  another 
argument.  For  simplicity  we  consider  only  the  homogenous  system 
( 3. 1)  dy/dt  = A(t)y 

where  A (t)  < C^,  p > 1 but  A and  its  derivatives  need  not  be  of  moderate 
size.  Assume  that  we  want  to  solve  this  system  by  a standard  difference 
method.  Then  we  discretize  (3.1)  and  use  point  values  A(tQ).  Therefore 
it  is  reasonable  to  demand,  that  the  equations 

dw/dt  = Aqw,  Aq  - A(tQ) 

with  constant  coefficients,  locally  describe  the  system  (3.1).  For  the 
numerical  technique  to  work,  it  is  essential,  that  the  solutions  of  (3.1) 
do  not  grow  too  fast.  Thus,  we  demand  that  there  are  constants  KQ,  c of 
moderate  size  such  that 


for  all  fixed  values  tQ.  (3.2a)  implies  that  the  eigenvalues  k of  AQ 
satisfy  the  condition 


Real  \ < c . 

Another  demand  is  that  the  solutions  do  not  oscillate  too  rapidly.  This 
leads  to  the  condition 

(3.2b)  |lm  k|  < p|  Real  \|  + c,  p,  c constants  of  moderate  size  . 
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We  shall  show  that  the  conditions  ( 1.  2a)  and  (3.2b)  naturally  lead 
to  matrices  which  are  essentially  negative  dominant.  We  have 


Theorem  3.1.  Let  3^(Kq,c,p)  denote  the  family  of  all  n Xn 
matrices  for  which  (3.  2a)  and  (3.  2b)  hold,  and  let  3(c,p,6)  denote 
the  family  of  n x n matrices  which  are  essentially  negative  dominant. 

For  any  fixed  KQ,  6 there  is  a universal  constant  such  that  for 

every  matrix  A < 3j(Kq,c,  p)  there  is  a nonsingular  transformation  S with 

IS  | |S-1!  < Kj  and  S_1AS  t 3(c,  p, 6)  . 

Proof.  In  [ 6 j we  have  shown  that  there  are  universal  constants 
*11,  ^12  suc^  t^lat  f°r  every  A t 3^  there  is  a nonsingular  transformation 


I T I I T-1 | < K 


T l(A  - cI)T 


b b 

11 

0 b 


I 

12 

b 

In 

_ b 

b 

22  23 

2n 

• 

0 

b 

nn 

X lb. . I < K | Real  b. . | Real  b..<0,  i = 1 2,  . . . , n . 

.tj  ij  12  li  n ~ ’ ’ 


llm  b. . I < p Real  b. , I + c,  i = 1,  2,  . . . , n 

li  ~ li  ’ ’ ’ 
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Therefore  we  can  choose  a diagonal  scaling 


h 


D = 


0 


0 


\ 


0 


d > 0 
i 


v° 


such  that  S ^AS,  S = TD  is  essentially  negative  dominant,  i.e.  belongs 
to  5^.  Here  j D i |D  only  depends  on  6,  K . 

Now  assume  that  A(t)  e 3^  for  every  fixed  t and  consider  the 
system  (3.1)  in  a neighbourhood  of  a point  t . Let  SQ  denote  the 
transformation  of  the  last  theorem.  Introducing  a new  variable  y = SQV 
we  obtain 


i 3.  3)  dy/dt  = B(t)y  , B(t)  = S^Aft^  , 

where  B(  tQ)  is  essentially  negative  dominant.  If  A(t  ) shall  represent 

A(t)  in  a whole  neighbourhood  of  tQ  then  B(t)  must  be  essentially 

negative  dominant  in  a neighbourhood  of  t^  provided  we  change  6 to  ~~  6 

and  c to  Zc  . Thus,  we  can  divide  our  interval  of  integration  into 

subintervals  t < t < t.  . and  in  every  subinterval  transform  A(t)  into 
l = = l +i 

an  essentially  negative  dominant  matrix  by  a transformation  S.  for  which 

Is  i |S.  * I is  uniformly  bounded.  A more  precise  description  is  given  in 
i i 

Theorem  3.2.  Consider  the  system  (3.1)  and  assume  that  there  are 
constants  K,  c,p  such  that  the  matrices  A(t)  « 3f^(K,  c,  p)  for  every 
fixed  t.  Assume  also  that  there  is  a constant  K such  that  for  all  fixed  t 

I A_1(t)dA/dt  | < K , A = A - ( c + 1)1  . 
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Then  one  can  divide  the  interval  of  integration  into  subintervals 
ti  < t < t.+1,  t.  +1  - t.  > n > 0,  r,  = const. , 

and  in  every  subinterval  there  is  a transformation  S.  with  |s.  I |s.H 

1 1 i 

uniformly  bounded  and  S.  ^A(t)S.  t g^(c  + 1 , 2 p ,6).  Furthermore  6 
can  be  choosen  arbitrarily. 

Proof.  Let  tQ  be  a fixed  point.  Without  restriction  we  can  assume 
that  A( tQ ) is  essentially  negative  dominant.  Otherwise  we  would  use 
the  transformation  of  theorem  3.1  to  obtain  the  system  (3.  3).  Proceeding 
as  in  lemma  2.  5 

A(t  - tQ)  = A(t0)E(t  - tQ) 

with 

E(0)  = I and  | d E/dt  I = |A_1dA/dt|  < 
and  the  theorem  follows. 

If  the  system  does  not  satisfy  the  above  conditions  then  using  a 
numerical  procedure  directly  can  be  very  dangerous.  However,  we  shall 
show  that  we  can  obtain  the  normal  form  by  combining  the  transformation 
of  the  dependent  variable  with  a local  stretching  of  the  independent  variable. 
We  divide  the  interval  of  integration  into  subintervals  t.  < t < t. 
and  construct  in  every  subinterval  a nonsingular  transformation  S.  and 
a stretching 

t - t.  = o.(t  - i),  a.  = t.  , - t.,  i < t < i + 1 , 

l i ’ l l+l  l 
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such  that  the  transformed  systems 
13.4) 


dv/dt  - a.S.  lA(t)S.v 


are  essentially  negative  dominant  and  the  functions  (1.8)  are  of  moderate 
size.  Then  the  solutions  of  (3.  4)  are  smooth  in  the  interior  of  every 
subinterval  i < t < i + 1,  but  boundary  layers  could  appear  at 
t 1,  2,  ...  . By  lemma  2.  6 we  can  avoid  these  by  demanding  that  the 
stretching  constants  u do  not  increase  too  fast  and  that  the  transforma- 
tions S do  not  change  too  fast,  i.e.  there  are  constants  d > 1 and 


Kq  > 1 of  moderate  size,  such  that 


3.5)  “i+A-d-  max(|s^S.+1l,  I S^S.  I ) < KQ  . 

We  shall  now  describe  the  details  and  start  with  the  scalar  equation 


,-l 


(3.  6) 


edy/dt  = -fy,  y(0)  = yQ,  t >0  , 


where  p is  a natural  number  and  e > 0 a small  constant.  Its  solution  is 
given  by  y = axp(tP+Vfc(P  + l))yQ,  it  is  a function  of  t/e^P 
In  a neighbourhood  of  t = 0 we  introduce  into  (3.6)  a new  variable  t 
by  t = at,  a = const.  > 0 and  obtain 

P+1  TP 


(3.7) 


dy/dt  = - 


— y = au(t)y  ■ 


We  determine  the  largest  a such  that 


(3.8) 


max  (min(l,  I a^  (t ) | ) I da^(t  )/dt  I < K . 

0 <t<l 
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* 


Here  K > 1 is  a threshhold  constant  of  moderate  size.  For  simplicity 
we  assume  that  K - p and  obtain 


pu 


P'l 


Ke, 


i . e.  u 


Thus  we  have  determined  the  stretched  variable 


l/(p+l) 

for  0 < t < e F = t, 


(3.8)  guarantees  that  the  solution  of  the  transformed  differential  equation  (3.7) 
is  smooth  for  0 < t < 1.  (In  general  one  has  to  include  more  derivatives 
in  the  expression  1 3.8).  However,  for  the  equation  (3.6)  the  inequality 
(3.8)  insures  that  a^(t)  and  all  its  derivatives  are  of  moderate  size. 

This  is  quite  common.) 

We  determine  now  the  stretched  variable  in  a neighbourhood  of  t - t^. 
Let  t = t + a(t  - 1),  then  the  differential  equation  becomes 

(t  + (t  - 1))P 

dy/dt  = “ y = a^(t)y  • 

The  obvious  modification  of  (3.8)  gives  us  a = t^  and  we  have  determined 

the  stretched  variables  for  0 < t < t = 2t^.  This  process  can  be 

continued.  After  n steps  we  have  obtained  the  stretched  variables  for 

0 < t < t = 2n_it  ■ It  is  clear,  that  the  interval  length  corresponds 
= = n 1 

exactly  to  the  behavior  of  the  solution  of  (3.  6). 

If  the  given  system  of  differential  equations  (3.1)  is  essentially 
negative  dominant,  then  we  can  use  the  corresponding  procedure.  In  the 
simplest  case  the  condition 


(3.9) 


max  (min(l,  I a^  (t ) I ) I da  (t  )/dt  I ) < K,  i,j  - 1,  2,  . . . , n , 
i <t  < i+ 1 


' 


! 
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determines  a stretching 


( 3. 10)  t - t - o (t  - t .).  t , - t 1 . 

l i l i+l  l 

Applying  (3.9)  to  (3.1)  gives  us 

( 3. 11)  dy/dt  a.A(t  )y  . 

If  u < 1 then  also  ( 3.11)  is  essentially  negative  dominant.  If  a.  > 1 

l — i 

then  this  need  not  be  so  and  we  have  to  limit  u to  a value  u.  > 1 

i i 

such  that  ( 3. 11)  satisfies  (1.  3a),  (1. 4a)  and  >1.  4b). 

If  the  system  (3.1)  is  not  essentially  negative  dominant,  then  we 
have  to  transform  it  to  that  form.  This  can  be  done  in  the  following  way. 
For  t - 0 we  use  the  Q - R method  to  construct  a unitary  matrix  U(0) 


t - pt,  0 < t < 1,  such  that  the  eigenvalues  of  B^O)  = |3B^(0)  satisfy 
these  conditions.  (This  procedure  is  not  very  satisfactory  and  shall  be 
improved  in  another  paper.) 


The  last  step  consists  of  a diagonal  scaling 


[ 


such  that 


(3.12)  DOlB2(0)D0  = B3(0) 

satisfies  the  condition  (1.4b). 

Unfortunately  i D I | D ^ | is  not  always  of  moderate  size.  In  that 
case  we  proceed  in  the  following  way.  Let  t > 1 be  a constant.  We 
shall  divide  the  eigenvalues  (3k.  of  B (0)  into  classes  N..  The  class 

Id.  J 

is  defined  by 

1)  (3k  . t if  I |3k . I £ c , 

2)  (3ku  e N1  , 

3)  (3k.  € N if  |3k.+^  « and  I k.  1/ 1 k^  I <t  . 

If  Nj  does  not  contain  all  eigenvalues  then  there  is  one  eigenvalue 

with  (3k  , ( N,  but  (3k  t N,.  Then  the  class  N_,  is  defined  by  the 

p+1  1 p i 2 

last  two  rules  with  (3k  replaced  by  (3k  • The  other  classes  are 

n p 

constructed  correspondingly. 

Now  we  construct  a transformation  H such  that 
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where  every  B is  of  the  same  form  as  B,(0)  and  its  eigenvalues  consist 
))  *■ 

of  the  class  N . Then  we  apply  a diagonal  scaling 

j 

fa  0 • 0 \ 

0 D 0 . . . 0 

D.  = 

1°  ■ ' 0 Drr  / 

such  that 

(3.13)  B3(0)  = D~1H~1B^(0)D1H 

is  essentially  negative  dominant.  If  I D H I I ( H ) *1  < I DQ  I I DQ  I we  use 

(3.13)  otherwise  we  return  to  (3.12). 

In  this  way  we  have  constructed  a transformation  v = SQ*y  and  a 
stretching  t = pt,  0 < t < 1,  such  that  the  system 


has  the  property  that 
there  is  an  interval 
if  we  replace  c by 


dv/dt  = [3S()1A(t)S0v 

pS~1A(0)SQ  is  essentially  negative  dominant.  Then 
0 < t < t in  which  the  same  is  true  for  pSQ  A(t)SQ 
2c  and  6 by  ~ 6.  Let 


t = 


t if  t > 1 

t/t{  if  ^ < 1 


then  the  system 

(3.14)  dv/dt  = a0S()1A(t)S0v 

is  essentially  negative  dominant  for  0 < t < 1. 


At  t 1 we  modify  this  process  somewhat.  By  (3.  5)  the  stretching 


shall  not  increase  too  fast.  Also,  in  many  cases  is  a good  approxima- 

tion for  Sj.  Therefore  we  perform  at  t = t = a preliminary  scaling 

*•  -1 

t - t - cic/g(t  - 1),  1 < t < 2,  and  the  transformation  u - y to  obtain 

(3.15)  du/dt  = da^S^ A( t)SQu  . 

If  (3.15)  is  essentially  negative  dominant  then  we  use  (3.15).  Otherwise 
we  use  the  process  as  described  earlier  to  construct  a matrix  T(t^) 
such  that 

(3.16)  (3da0T"1(t1)SQ1A(t1)S0T(t1) 

is  essentially  negative  dominant.  If  j T I | T is  of  moderate  size  then 
we  proceed  as  earlier.  Otherwise  we  return  to  (3.14)  and  replace 
and  t by  aQ/2  and  t^/2  respectively  and  restart  the  process. 
Eventually  the  procedure  will  stop  because  in  the  worst  situation 
d(a0/2P)S0^A( t ( 2p)SQ  will  be  of  moderate  size.  Then  we  can  choose  T = I, 
i.e.  SQ  = Sj  and  proceed  as  earlier.  In  this  way  we  determine  t 
and  the  following  t^  and  construct  a system  which  is  piecewise  negative 
dominant.  The  new  system  is  treated  as  stated  earlier  (see  (3.9)). 

We  shall  now  illustrate  the  technique  for  a simple  example.  Consider 
the  differential  equation 

ey"  + (a(t)y)'  + b(t)y  =0,  t > 0 

with  initial  conditions 

y(o)  = y0,  y'(0)  = . 
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^ 


4 


Here  e > 0 is  a small  constant  and  a,  b with  a > 0,  b < < 0 are 

smooth  functions  of  moderate  size.  We  write  the  above  equation  as  a 
first  order  system  by  introducing  new  variables  y^  = y,  dy^  ^/dt  = b(t)y 


and  obtain 


(3.17) 


dy  d 

l^\ 

’ 

dt  ~ dt 

[ 

\ 

a(tl 

e 

b(t) 


|\  / / 1 \ \ 


(1) 


(2) 

/ 


= A(t)y 


If  aft)  > 1 and  ib(t)  I < c then  the  system  is  essentially  negative 
dominant  and  nothing  needs  to  be  done.  Otherwise  we  transform  the 
matrix  A for  t = 0 to  upper  triangular  form  employing  a unitary  transforma- 
tion U - U(0) . The  eigenvalues  \ of  A(0)  are  the  solutions  of  the 
characteristic  equation 

e\2  + a(0)\  + b(0)  = 0 , 


i . e . 


x a(0) 

\ = - 2e  Sign  a ' 


kill)  , a LQ.) 

e 2 ’ 

4e 


a(0) 

X.  , = - + sign  a • 

2 2e 


b( 0)  a (0) 
£ 4e2 


Observe  that  for  la  I » Ve  the  two  eigenvalues  have  the  form 

a{0j  v biOj  biOj 

K ~ “ , A.  = " 


a(0) 


Furthermore  b < 0 implies  that  \ , y are  real  and 

lyo)  I > |b(0)/e  I1/2  • 
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The  eigenvector  corresponding  to  V is  given  by 


e , = x , t = (1  + (b(0)/K  (0))  ) 

\b(0)A.(0)/ 


and  the  desired  unitary  transformation  has  the  form 


mi 


U(0)  = T 


. b(0) 


Bj(t)  - U (0)A(t)U(0) 


K ( 0)g(t)  - - 

1 r 


dit)  b(  0 ) 
e^(0)  e\I(0)< 


with  d(t)  = b(0 )(a(t)  - a{ 0)),  and  g(t)  = 1 + 


e\(0) 


Thus  d(0)  = 0,  g( 0)  = 1 


and  B (t)  is  upper  triangular  for  t = 0.  The  condition  (1.  3a)  is  always 

satisfied.  However,  (1.4a)  does  not  need  to  hold,  because  \^  > 0 
and  \ = 0(lAe)  when  a(t)  = 0.  To  satisfy  condition  (1.4a)  we 

A 

perform  a stretching  t = pt,  where 


i lf  -M2L  < i c 

“ t^lO)  - 2 • 


ctV0)  _bi»L  i 

2b(0)  1 t\  (0)  > 2 c ' 
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I 


Then 


B^t)  = pBx(t)  = 


^KA(0)g(pt) 


V 


The  next  step  is  to  apply  a diagonal  scaling 


-fi  \ 


)3d(pt)  bb(0) 


eV(O)  eK  (0) 


I 


I 2 


D(0) 


z\{0) 


1 J 


Let  v = D y.  Then  (3.17)  becomes 

/ 


3. 19) 


dv 

dt 


p\1(0)g(pt) 


2 2 v 

^ e ^(0) 


2 V0) 


23 d(Bt ) i3b(  0 ) 


e^(0) 


v . 


Denote  by 


(3.  20) 


f ceX.  (0) 


* ceK^O) 

*1  - 2pb(0)  = ' 


2b(0) 


if 


if 


b(0)  „ 1 

< — c 


e^(0)  - 2 


^>ic 

e^(0)  2 ’ 


A A A 

the  largest  value  such  that  for  all  t with  0 < t < t 


1 


(3.21) 


a(0)  - a(  Bt  ) > X 
e\(0)  ~ “ 4 


; 2gd(pt) 
1 e ^(0) 


^Ic’ 


A A AW 


and  introduce  a new  variable  t by  t = t^t.  Then  (3.19)  becomes 
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w 


I 


and  the  last  system  is  essentially  negative  dominant  for  0 < t < 1.  The 
conditions  (3.20)  and  (3.21)  are  satisfied  if  we  choose 


(3.23)  tj  = Ptj  = CjeUUO)  | . 

Here  c^  is  a constant  which  depends  on  c,  b and  da/dt  but  not  on 
e.  Also,  the  conditions  (3.9)  are  essentially  satisfied  and  no  further  stretching 

is  necessary.  Thus,  the  interval  length  0 < t < t^  is  given  by 

t - const,  e |\j(0)  |.  In  general  we  get 

t.  - t.  = const,  e I \(t.)  I . 
l+l  i l 

There  are  two  different  situations 

1)  a(t)  > aQ  > 0 everywhere.  Then  e I = I a/b  I and  the  interval 
length  never  becomes  very  small. 

2)  a(t)  = 0 for  some  t = t . Then  e U(t ) I = 0 (Ve)  and  the 
interval  length  deminishes  exponentially  to  0('v/e)  when  t approaches  t. 

This  is  completely  in  accordance  with  the  behavior  of  the  solution  of  the 
differential  equation. 
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4.  Difference  approximations. 

In  this  section  we  consider  the  system  (3.1) 

(4.  1)  dy/dt  = A(t)y  f F(t),  t > 0 , 

and  approximate  it  by  general  difference  approximations.  We  start  with 
multistep  methods, 


r r 


(4.  2) 

L[  u]  = ^ (a  I + 

j = -l 

II 

1 

D 

1 

4-> 

ro 

-k  l P,F(tv-j>  ■ 

j = -l 

where 

a = 1 and  p ^ < 0 

. We  write  (4.2) 

as  a one-step  method 

14.  3) 

v - B(t 

v +1 

, k)v  + kF  , v 

V V V 

= 0,1,2,  • • • 

with 

ft'  \ 

I Uv+\ 

, F^U+k^Alt^lf1 

(-  A 

V - 

V 

= 

0 

V 0 i 

B(t,k)  = 


f E0(t,k)  ....  Er(t,k)^ 

I 0 . . • o 

0 10..  0 

^ 0 ..01  0 ^ 


where 


E (t,  k)  = “(I  + k(3  A(t  + k))  V.I  + kp.A(t  - jk)),  j = 0,1,  . . .,r 
j ”1  J J 
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The  aim  of  this  section  is  to  prove  theorem  1.3.  For  this  we  need 


the  following  two  well  known  lemmata. 

Lemma  4.1.  Consider  a class  X of  problems  and  assume  that  the 
approximation  is  stable  for  this  class.  Let  A(t)  be  the  defining  matrices 
and  consider  the  class  X of  problems  where  A(t)  is  replaced  by 
Alt)  = A( t)  + B(l).  If  the  B(t)  are  uniformly  bounded  then  the  approximation 
is  also  stable  for  the  class  X • 

Lemma  4.2.  Consider  a class  X of  problems  and  assume  that 
there  are  constants  r|  > 0,  K > 0 such  that  for  every  problem  in  X : 

1)  the  interval  of  integration  can  be  divided  into  subintervals 

t.  < t < t.  , . with  t.  - t.  > n 
l ~ ~ l+l  i+1  i ~ 1 

2)  there  is  a transformation  T.(t)  in  each  subinterval  which  satisfies 

It.  | It"1  | + |dT./dt|  < k 

3)  the  approximation  (4.2)  is  stable  in  each  subinterval  for  the 
transformed  problem 

dy/dt  - Ay  A = T_1AT.  + T_1dT  /dt  . 

l l i i 

If  the  local  stability  constants  Kj,  Kg  are  uniformly  bounded  then 
the  approximation  is  also  stable  for  the  class  X • 

We  make 

Assumption  4.1.  We  consider  a class  ^(fi)  = X(n,  c,  p,  6,  K)  (see 
definition  1.  2)  for  which  c = 0 and  1-6  is  sufficiently  small.  Also, 
the  approximation  (4.2)  satisfies  assumption  1.2. 
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This  assumption  is  no  restriction.  1)  If  c * 0 then  every 

Ai  ?r(S2)  can  be  written  in  the  form  A = where  the  are 

uniformly  bounded.  By  lemma  4.1  we  can  neglect  B . 2)  Let  A( t^) 

satisfy  the  inequalities  (1.  3a),  (1.  4a)  and  (1. 4b)  with  c = 0.  By  lemma  2. 2 
A(t  )t 

le  111  i.e.,  the  relation  (3.2a)  is  satisfied  with  K0=1»c  = 0. 

Furthermore,  for  every  eigenvalue  k of  A(tg),  there  is  an  a. . such  that 

n 

U - a..  I < X la. . I 1 (1  - 6)  I Real  a I . 

11  ~j  = i 

j*i 

Therefore, 

I Real  X.  I > 6 |Real  a I,  |lm  \ | £ I Im  a. . I + (1  - 6)  I Real  a..  I £ 

< ( p + 1 - 6)  | Real  a I £ 6_1(  p + 1 - 6)1  Real  \ I , 

and  the  relation  (3.2b)  holds  with  c = 0 if  we  replace  p by 
p - (p  + 1 - 6)/6.  Thus,  if  A(t)  e 7r{U,  0,  P,  6,  K)  then  A(tQ)  f ^(1,0,^). 
Furthermore,  A(t)  - I is  negative  dominant  and,  by  lemma  2.4,  the 
logarithmic  derivative  (A  - I)  dA/dt  is  uniformly  bounded.  Therefore, 
by  theorem  3.2,  we  can  transform  the  class  ft'(n)  to  another  class  ’r(n) 
where  1-6  is  as  small  as  we  like.  The  transformation  is  piecewise 
constant  for  every  A and  |t.  I |t.  1 | is  uniformly  bounded.  By  lemma  4.2 
we  need  only  consider  this  new  class. 

To  simplify  the  notation  we  shall  use 

Definition  4. 1.  A(t)  is  weakly  negative  dominant  if  (1.  3a),  (1.4a) 
and  (1.  4b)  hold  with  c = 0. 


We  shall  also  simplify  the  difference  approximation.  For  c - 0 


the  matrices  Aft)  - 1,  A(t)  < 7r(ii)  are  negative  dominant  and,  by  lemma  2.5, 
' 4.  4)  A(t  - jk)  - 1 - (Aft)  - I)(I  + jkC  (t,  k))  , 

where  the  C are  uniformly  bounded.  The  matrices 

1 

4.  5)  E (t,  k)  - (I  + kp_1A(t))“1(ad  + kp.A(t)) 

are  uniformly  bounded  and  have  uniformly  bounded  derivatives.  Further- 
more, by  ( 4.  4), 

E (t,  k)  = E (t,  k)  + kC.(t,  k) 


where  the  C.(t,  k)  are  uniformly  bounded.  Therefore,  by  lemma  4.1, 
we  can  replace  (4.  3)  by  the  difference  approximation 
(4.6)  v - B(t  . k)v  + kF  , v = 0,1,  2,  ... 

V +1  V V V 

where  B has  the  same  form  as  B with  E replaced  by  E . . 

We  consider  now  some  special  cases.  Let  k be  fixed  and  denote 
by  Tq  , the  class  of  problems  where  A is  negative  dominant  and 
! kA  i < d for  all  t.  Here  d is  a sufficiently  small  constant.  By  (4.  5) 

E.  = a. I + p.kA  + kAQ.kA,  + p.  , 

where  0 are  analytic  functions  of  kA.  Therefore  we  can  write  the 
3 

matrix  B(t,  k)  in  the  form 


B = Bq  + ABL  + AQA 


where 


| 
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/ k A(  t ) 0 . ...  0 \ 

0 kA(t)  0 ...  0 


A - A(t)  - 


tion  (4.2)  is  stable  in  the  usual  sense.  Therefore,  there  is  a nonsingular 
transformation 


ftni 


T = 


W \ 


such  that 


Vr+ll1  ‘ ' ' V+lr+l  / 


lD  1 


tbqt 


-1 


o d2  0 


0 \ 

0 


V 0 


0 D / 
s 


where 


and 


n 


- k.(0)I,  |k.(0)  I = 1,  j = 1,  2,  • • s - 1,  k * k.  for  j * i , 


J J 


Id  I <1-0-,  cr  = const.  > o . 
S 


\ 0 

0 

kAlt )/ 

5 

fV  “i1 

. . . 

“rM 

(V 

K1  ••• 

Bo^ 

1 0 

0 

’ V 

0 

• 

. 

0 

1 0 

0 I 

0 ) 

l 0 

. 

0 ) 

is  an  analytic 

function  of 

kA.  By  assumption  1.  2, 

the  approxima- 
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We  can  consider  the  matrix  AB^  as  a perturbation  of  BQ.  There- 
fore, we  can  use  the  following  theorem  of  B.  Engquist  [4|. 

Theorem  -1,1.  There  is  a transformation 

H T + AR 

where  R is  an  analytic  function  of  kA  such  that 
4.  ')  HBH"1  = D + AQA  . 

Also,  Q is  an  analytic  function  of  kA  and 


(«£  0)(I  + kg^A) 


k,(0)(I  + kgzA)  0 


. . . 0 . 


D + kg  A 
s s 


We  can  now  prove 

Lemma  4.  3.  Consider  the  class  of  problems  7) ' ^ For  sufficiently 

small  d and  k there  are  constants  n > 0,  > 0 such  that  the 


solutions  of  (4.  6)  satisfy  the  estimate 

qq(t  ) 

(4.9)  lv(tv)|  £ Kls(e  V |vn 


/ _ 

y | v I + max  I a"  (t.)F(t  ) 
0<j<v  J J 


(4.10) 


q t ) = ),  a t,  k,  a(t  ) = max  Real  a,,  t.)  . 

v n L-J  , 1 j . 11  j 

0<  <v-l  J ' i 


Proof.  Let  w = H v . then  the  equation  (4.6)  can  be  written  as 

■»" V V V 


(4.11) 


H H”  w , = (D  + A 0 A )w  UH  F . 

V v-fl  v+l  V v V V V V V 


By  lemma  Z.  5 


kAft  + k)  = kA( t)(I  + kC(t,  k)) 


where  C is  uniformly  bounded.  Therefore 
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H H 


i +kA  P 

i'  1/ 


where  P is  also  uniformly  bounded,  f urthermore, 
v 

kH  F A G , G = kH  A~V  , 

V V V V V V V V 


with 


(4.12)  IG  I < const. 


r -1  r.  -1 

V l A ( t )F(t  ,)l  < const.  X i A~  ( t ,)F(t 
— 1 , v'  v-y  - . v-)  v-) 

j = -l  J = -l 


We  first  consider  the  homogenous  equation  (4.11),  i.e.,  F = 0.  For 
sufficiently  small  d,  the  row  sums  R ^ of  the  right  side  of  (4.11) 
can  be  estimated  by 


4.13) 


R < \ - - cr,  for  the  rows  corresponding  to  D , 
12  s 

R < lick  Real  a + g k Im  a . . I + ( 1 - 6 )g  , k i Real  a . 

1 - j ii  j n ) u 

+ const,  d ( k Real  a..  + k Im  a. . I < 
a n 


< 1 - |g  k Real  a I + const.  d(l  + 6)  |k  Real  a. . i £ 

— j ii  ii 

1 - — g k I Real  a. . I otherwise  . 

2 ] li 

For  the  row  sums  Rz  of  the  left  hand  side  we  obtain,  correspondingly 

R > 1 - const,  k^  i Real  a. . i 


and  the  lemma  follows  from  F = 0. 

Now  assume  that  wQ  = = 0.  Consider  (4.11)  for  0 < t . < t^ 

and  let  w^1  ^ = max  I w^  | = ||w  || . Using  row  a of  (4. 11)  with  v + 1 = r 

* U J 

we  obtain  from  (4.13) 


1) 


The  row  sums  of  a matrix  (a. .)  are  defined  by  I a . I , i 1,  2,  . . . , n. 

1J  J = 1 1 
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or 


0(k))||w||  < const,  k ! Real  a. . I II G |i  , 
- n 


g.(l  f 0(  k))k  | Real  a w ii  < const,  k I Real  a I |j  G li  , 

which  proves  the  lemma. 

The  next  special  case  is  treated  in 

Lemma  -1.4.  Let  cl^,  be  constants  with  0 < d^  < d^,  and  k 
be  iixed.  Consider  all  problems  belonging  to  ?r(Q)  for  which  the  eigen- 
values of  the  matrices  kA  belong  to  C> , (see  definition  1 . 2) . Then, 

dr  a2 

there  are  constants  and  a,  with  0 < u < 1,  such  that  the 

solutions  of  (4.  6)  satisfy  the  estimate 


4.14) 


i v(t  ) i < K 
v 2S 


(a  + 0(k))  i vq  ' + max 
0 <)  < v 


A~1(t.)F(t.) 


1 - <r  + 0(k) 

Proof.  The  eigenvalues  k of  B have,  by  assumption,  the  property 


that 


< 1 - T . 


Furthermore,  B and  dB/dt  are  uniformly  bounded.  Therefore,  we  can 
without  restriction  assume  that  j B | < 1 - j t.  Otherwise,  we  can  find 
a transformation  T with  I T i | T + IdT/dtl  uniformly  bounded  such 
that  TBT  * has  this  property.  Also, 


r i r 

-L.  .i  \ 


.-1, 


k J F I < const.  |A  (t  )F(t  ,)j  < const.  X iA  *(*  .)F(t  .)  I . 

V V v-l  - , v-l  V — 1 

J = -l  J = -l 


Therefore,  (4.14)  follows  from  known  principles. 


We  can  now  prove  theoiem  1.3  for  the  case  that  the  approximation 
(4.2)  is  strongly  stable  for  the  class  T(c2)  of  scalar  equations . By 
lemma  4.  2 we  need  only  prove  local  stability.  We  consider  the  homogenous 
system  4.1  in  the  neighbourhood  of  a fixed  point  t^ . By  lemma  4.1 
we  can  assume  that  A t)  is  negative  dominant.  Let  k be  fixed  and 


d > 0 a constant.  Then  we  can  write  A(  t^ ) in  the  form 


AUq) 


^ A (t  ) A ,(t)'' 
110  120 


W A22U0 


I a 


i la 


nn 


where  A (t  ) is  the  largest  submatrix  for  which  2|kA  (tjl  < d.  For 

every  eigenvalue  \ of  A^  there  is,  by  the  Gershgonn  estimates,  a 

diagonal  element  a with 
u 

n 

I X.  - a. . I < ^ la  I < (1  - 6)  | Real  a.  I,  i.e.,  I \ | > 6 | a I . 
n — j ij  _ n ii 

j*i 

Furthermore,  the  row  sums  r of  A satisfy  the  estimate  kr  > ~ d. 

ill  l 2 

Therefore 


Ika. . I = kr.  - k ' Ja..j>~’d-(l-6)jka..|. 

n l ij  _ 2 n 

Thus  lka..|  >"~d  and  the  eigenvalues  \(A  ) of  A satisfy  the  estimate 

ii  4 1111 


lk\(A  ) | >~db. 

1 1 4 

By  lemma  2.  5,  there  is  a neighbourhood  It  - t I < r|  > 0 a constant 
independent  of  A and  tQ,  such  that  the  matrices  A..(t),  i = 1,2 
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have  the  property 


IkA  (t)|<d,  IkMA,  , )(t)  I > d.  d >~do. 

C.L  11  1 1 r> 


Now  consider  the  approximations  (4.6) 
(4.  l 6) 


v[\\  = 1 l(t  , k)vl  11  + kFl 

-f-  1 


V V 


.4. 17) 

for  the  subsystems 


v^  2*  = 2^(t  , k)v^  ^ + kF^2 

V +1  V V 


(4.  18) 
(4. 19) 


dy'  1 Vdt  = Ai  ( y*  1 ^ + F^  1 ^ 


dy^/dt  = A2^y^2*  + F^ 


respectively.  We  choose  d so  small  that  (4.17)  satisfies  the  estimate 
(4.9).  If  6 --  1,  then  the  eigenvalues  of  A converge  to  the  eigen- 
values i of  A with  J k\  | > d • By  assumption  1 . 2 these  eigenvalues 
belong  to  Q . Therefore,  the  eigenvalues  k of  B^(t,  k)  satisfy 

dr 

the  inequality 


Ik  I < 1 " r for  I t ~ t I < g 

provided  we  have  choosen  l - 6 sufficiently  small.  Thus,  the  solutions 
of  the  system  (4.  1 6)  satisfy  the  estimate  (4.  14). 

We  can  now  write  the  differential  equation  (4.1)  in  the  form  (4.  18) 
and  (4 . 1 9)  with  F^  ^ = A^  ^ 2 ^ , F^  2 ^ = A y^  ^ . The  difference 
approximation  (4. 6)  can  then  be  represented  by  (4. 16),  (4.  17)  and,  for 
sufficiently  small  k,  the  estimates  (4.14)  and  (4.9)  give  us 
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||v^  1 hi  £ fv^1  ^ I + const.  I|A[  Ja^  ^ II  llv^  II  , 
||v^  2 hi  £ Iv^2  ^ | + const.  IIa,~|a2i  ||  ll  1 hi  , 


max  |u 
0 £ j £ v J 


Therefore,  the  approximation  is  stable  provided 


is  sufficiently  small.  This  proves  theorem  1.3  when  the  approximation 
is  strongly  stable. 

The  proof  of  theorem  1.3  for  the  case  that  the  approximation  is  only 

stable  is  more  complicated.  We  start  again  with  a special  case.  Let 

k be  fixed  and  let  Vr  , denote  the  class  of  problems  for  which  A is 

d,  00 

negative  dominant,  A 1 is  weakly  negative  dominant  ( see  def.  -4.1), 
and  |(kA)  1 | £ d.  Then  we  have 

Lemma  4.  5.  Consider  the  class  of  problems  % , ..  • For  sufficiently 
d,  * 

small  1 - 6,  d and  k there  is  a constant  K > 0 such  that  the 
solutions  of  (4.  6)  satisfy  the  estimate 

(4.  20)  |v(t  ) I < K ((1  + kdt  ) |v  | + k(l  + dt  ) max  I F(  t . ) |)  . 

v — 3S  vO  v r\  s'  s 1 

0 £ ] £ V 

Proof.  By  (4.5),  the  E^(t,  k)  are  now  analytic  functions  of  (kA) 
and  we  can  write  the  matrix  B(t,  k)  in  the  form 

B(t,  k)  - Bq  + AB  + AQA 

where  0 is  an  analytic  function  of  (kA)  1 , 


(kA)’1  0 ...  0 

0 (kA)‘‘  0 ...  0 

0 . 0 . . . (kA) 

and  BQ,  are  constant  matrices  of  the  same  form  as  BQ,  B . By 
assumption  the  approximation  is  stable  in  a neighbourhood  of  \k  - . 

Therefore,  we  can  apply  the  same  process  as  used  in  lemma  4.  3.  Using 
theorem  4.  I again  we  obtain  an  equation  of  type  (4.  1 1 ), 

II  H w = (D  + A Q A )w  +kHF 

V V + l U+l  V V V V V V 

where 

w =Hv  H = T + AR,  R is  an  analytic  function  of  (kA)  1 , 

v V v 1 

and 

k (0)(I  + g j (kA)’ 1 ) ....  0 ...  0 

0 K^(0)(I+g2(kA)_1)  0 ...  0 

0 ....  D + g (kA)’ 

s s 

By  assumption,  A is  negative  dominant.  Thus,  by  lemma  2.5, 
(kA(t  + k))'1  = (I  + kC)(kA(  t))~ 1 , 

and  therefore 

H H 1 = I + kP  A , with  P uniformly  bounded, 

v v+l  v v+l’ 

Without  restriction  we  can  assume  that  D only  consists  of  blocks 

k.(0)(I  + g.(kA)’1)  because  the  block  D +g  (kA)  1 has  no  influence 
j J s s 
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on  the  form  of  the  estimate.  Also,  we  may  assume  that  all  k^(0) 


= 1 


If  not,  we  can  introduce  a new  variable  w by 


/«: 


w = 
v 


0 \ 

0 


VO 


Then  we  can  write  (4.  20)  in  the  form 


w 


v / 

s- 


(4.21)  (It  kP  A )w  = (I  + GA  + A Q A )w  + kH  F , 

' ' V V + 1 V + 1 v V V V V vv 


where 


( V 0 


G = 


g2I  0 ... 


0 \ 

0 


V 


0 gs-fV 


We  want  to  estimate  w = A w . Multiplying  (4.2l)  by  A gives  us 

1/  y V v 

(4.22)  (I  + kAvPlv)wy+1  = (I  + GAv+A^Qv)wv+kAvHvFv,  w0  = %W0 

Equation  (4.  22)  is  of  the  same  form  as  (4. 1 l).  (Observe  that  A^Q^  = Qv  ^ 

Therefore,  the  estimate  (4.9)  holds  and  we  obtain 

|w  I < K ( |w.  I + k max  Ih.F,  I)  . 
v ~ 13  ° 0 < j < v 3 3 


We  can  now  write  equation  (4.  2l)  as 


2~ 


(4  23)  w = (I  + GA  + A Q A )w  + kH  F + k PvWv  + l ' 

v+l  ' V V V V V V V v V + i 
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> i 


quation  4.2i)  is  again  of  the  form  ( 4 . 1 1 ) . We  split  its  solution 


(2| 


into  two  parts  w w'  ' 1 + wl  1 where 

V V V 


4.  24)  w^  1 ^ (1  + GA  + A Q A )w^  1 ^ + k^P  w , w|. 1 ^ - w , 

V+l  v v v v v K v +1  ’ 0 0 

4.2  5 w^  il  t GA  lA  Q A )w^^kA  H A 'f  , wj/  ^ = 0 . 

l>  fl  l’  I'WV’l’  V V V V 0 

The  approximation  (4.23)  is  stable.  Therefore  we  can  estimate  its 
solution  by 


i w 


1] 


< lw.  j + kt  max  |P.w.|  < 


- ' 0 


0 < j < V 


1 1 


il  + const,  kdt  ) lw  1 + const,  kdt  max  If.  i . 

v 0 v n „ j 


0 <j£  v 


(4. 9)  gives  us,  tor  (4. 25), 


i w 


[2] 


< const,  k max  i A 1 F . i < const,  k max  | F 


j J 


0<j<v 


0 < j < v 

and  the  estimate  (4.20)  follows. 

We  now  consider  the  full  system  4.1.  Corresponding  to  our  earlier 

procedure,  we  can  write  the  matrix  A in  the  form 

f A A A \ 

11  12  13 


A = 


A2l  A22  A23 

\A3l  A32  A33  j 


, lanl>la22l> 


> a 


nn 


in  a neighbourhood  1 1 — t I < ig,  rj  > 0 of  every  point  t^,  where  A3 3* 


11 


are  the  largest  submatrices  with  IkA^I  £d^  and  l(kA 


V“1 


11 


< d.  . 


Then  there  are  constants  d^,  ^22  suc^ 


-58- 


(4.26)  d < k |\.(A22)  I < k Ia^!  < d22,  MA^)  eigenvalues  of  A^  . 


(4.27)  I A^A^.  I < const.  (1-6),  | kA  ^ |<  j kA  ^A^A  | < const,  d^i-6),  j-1,3 

(4.  28)  i A , A I < const.  (1  - 6),  |kA  I < const.  d22(l  - 6),  j = 1 , 2 . 

Without  restriction  we  can  also  assume  that 

IkA  . I £ const.  (1-6),  j - 2,  3 . 

Otherwise,  we  apply  the  transformation 

f1  S.Z  S1 3) 


0 0 I 


S = 0 I 0 , S - -A~  A , j = 2,  3 , 

’ lj  11  1 j’  ’ ’ 


and  a simple  calculation  using  (4.27)  and  (4.28)  shows  that  S 1 AS 
has  the  desired  property. 


Now  consider  the  approximation  (4.  6) 


(4.28) 


v“'  = Bli| 

V+l 


(t  ,h)v^'  + kF^1',  i = 1,2,3, 

u \J 


for  the  subsystems 


(4.29) 


dyf  l]/dt  = A..y[l],  i = 1 , 2,  3 , 


and  assume  that  the  estimates  (4.  20),  (4.  1 4)  and  (4. 9)  are  valid.  The 
full  system  can  be  written  in  the  form  (4.  29)  with 

= y . 

A ■' 

Then  stability  follows  in  the  same  way  as  earlier  from  the  above  estimates 
provided  1-6  is  sufficiently  small. 
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This  proves  theorem  | . j for  the  case  that  A~j  is  weakly  negative 


dominant . 


We  shall  now  show  that  one  can  always  transform  A]  ^ so  that  A^  J 
has  the  above  property.  There  is  one  case  for  which  this  is  immediately 


Lemma  T.t>.  Let  A be  negative  dominant  and  assume  that 
1/T  < la..  1/ 1 a I < t,  i * j,  i,  j = 1 , 2,  . . . , n , 


where  t > i is  some  constant.  If 


or  . 1 - 6 /,  2 

T I ~ < 1,  a = —7—  VI  + p 


Then  A is  weakly  negative  dominant  with  constants 


Proof.  Let 


p + o'  , a 

p < i , 6 < 1 - t 

K1  - 1 + 0-  ’ 1-  1—  (X 


'*11  0 


0 a22  0 ...  0 


0 a 


A_1  = (I  + A 1 ( A - A ))~ 1 A 1 = (I  + C)A_1  = D = (d  ) 


The  norm  of  C = (c. .)  satisfies 


lcl<^ 


Therefore,  we  obtain 


- Real  d > - Real  ~ - — : — g r > - (Real  )(1  - a) 
a _ a.  6 la.  I ~ a. . 

n n 11 

I Im  d I < I Im  ~ I + ~1 T < I Real  ~ I ( p + <r) 

li  ~ a.  6 a. . I — a.. 

u n li 

for  the  diaaonal  elements  d . of  D.  Furthermore,  by  (4.  30), 

ii 

n n n , 

d 'h  ic../a  I < ( t/  |a,.  I ) \ I c . I <TfT|Real  I < ~ ~|Reald.. 
id  ‘j  itt  » ” - u jf,  ~ au  ~1-’ 


This  proves  the  lemma. 

If  the  diagonal  elements  of  A ( are  all  of  the  same  order  of  magnitude, 
i . e . , t is  bounded,  then  we  can  assume  that  j is  weakly  negative 
dominant,  because,  as  we  have  seen,  we  can  make  1-6  as  small  as 
we  like. 

We  now  consider  the  case  that  A^  ^ contains  eigenvalues  of  different 
orders  of  magnitude  and  we  want  to  transform  A^  into  a block  form 
which  separates  these  eigenvalues.  We  need 

Lemma  4.  7 . Let 


A = ’ 


’ ail~  a22  ~"'~  ann 


be  a weakly  negative  dominant  matrix  and  assume  that  there  is  a constant 

t > i such  that  la  l/la  . , I > T 

rr  r+1  r+1 


for  some  r.  If  T(1  - 6)/(t  - 1)  is  sufficiently  small  then  there  is  a 


-fel  - 


transformation  S with 


max{ IS | , |S  | ] < 1 + T 


T.ll.  ~ f>)  K 

t - 1 K1 


such  that 


.4.  31)  A - S AS 


A11 


0 A, 


air  % hi 


the  matrices 


A.  - K — r 

11  2 T - 1 


are  also  weakly  negative  dominant  and  are  universal  constants. 

Furthermore,  the  elements  of  S are  analytic  functions  of 
a a a 


(4.  32) 


a a a 

E3 . . . p£ 

a - a a a - a . 
n jj  pp  li  jj 


, i = 1 ,2,  . . . , r;  j = r + 1 , . . . , n;  p > i,  p * q 


A = A+A-A  = 


A 

in  the 

o’ 

3 

/A 

1 

0 

+ (1  - 6) 

B 

1 1 

Bu' 

,o 

Azl 

iB21 

B22  J 

A 1 ° 


0 A- 


B 0 

1 1 


0 B_ 


w 


where 


an  0 


0 a,,  0 


r+1  r+1 

ar+2  r + 2 


For  S we  make  an  ansatz 


0 a 


S = I + (1  - 6 


0 T 


T2i  0 


Then  the  relation  (4.  31 ) is  equivalent  to: 

Bn+(1-6>  B12TZ!  = *„'  822  + (‘-6)  B2iT12  = ~B22' 

hlB,2  + T,  2 - Ai‘T12A2  = » ' 6>‘A;lT,2'B22  ' Ai‘B,lTl2»  > 
B21A71  - T2,  + = (1  - 6)<T2,51,A7‘  - B22T2,A7l)  ' 

Neglecting  terms  of  order  (1-6)  in  the  last  two  equations  gives  us  a 
linear  system 


A7‘B12  + T12  - A7‘T1ZA2  = ° 

Ba,A7‘  - T2!  + A2T21A7‘  = ° 


which  has  a unique  solution.  The  elements  of  T arc  of  the  form  (4.  32) 

U 

and  its  row  sums  are  bounded  by  V(T  - 1).  Therefore  the  lemma  follows 
easily  by  iteration . 

We  con  now  complete  the  proof  of  theorem  1 - 3.  Consider  the  matrix 


fa 


An(t) 


1 1 


Varl 


aiA 


rr  ‘ 


for  a fixed  point  t - t„.  Assume  that  the  a.,  are  ordered  such  that 
^ 0 u 

la22(t0)l  larrV!  ' 

1 /r 

Let  t | 7 >i  be  a constant.  We  divide  the  a,,  into  classes  7?.  by 
the  following  procedure. 


1)  au  < V 


l/r 


2)  a <77  if  a.  <7-  and  la.  , . ,/a..  <t 

jj  1 j-1  J-l  1 J-l  3-1  )j  1 

3)  If  V does  not  contain  all  eigenvalues  then  there  is  a last 


a < 77  with  a i V . 7?  is  then  defined  by  the  first  two  rules 

p-1  p- 1 1 PP  1 2 

with  a replaced  by  a 

1 1 PP 

This  division  of  the  diagonal  elements  of  A into  classes  at 

t - tg  defines  a division  into  classes  in  a whole  neighbourhood 

It  - t I < n of  to-  t-e--  a (t)  < ??,  if  a..(tQ)  < 7?..  By  assumption 

la”*da  /dt  I are  uniformly  bounded.  If  we  choose  n sufficiently  small 
33  jj 


then  the  a 

3) 


belonging  to  the  same  class  satisfy  the  condition  of  lemma  4.  6. 


Those  Leh  nging  ’ different  lasses  satisfy  the  condition  of  lemma  4.7 


provi  ie  : l - is  suffi  lentiy  small.  Repeated  use  of  lemma  4.7  shows 

that  we  in  !:  in  : A ^ to  blockdiagonal  form  where  every  block 

ati;  : • . ns  f lemma  4.6.  Furthermore,  the  transformation  S 

is  such  the-  S S 1 + iS/dt  is  uniformly  bounded.  Therefore 

theorem  1 5 follows  from  lemmata  4.  l and  4.2  and  the  earlier  stability 

results  - 

No  new  difficulties  arise  if  we  consider  approximations  of  type  (1.29). 
We  can  again  split  the  matrix  A into  three  parts  and  it  is  not  difficult 
to  prove  analogs  of  the  lemmata  4.3-4.  5.  Therefore,  theorem  1 . 3 als 
holds  for  approximations  of  type  (1.29). 
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